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ABSTRACT 


Viscous  drag  in  the  transonic  regime  over  an  axi-symmetric  body  with  a 
unique  aft  contour  surface  is  investigated.  The  forebody  is  composed  of  an 
arbitrary  ellipsoid.  The  unique  aft  contour  surface  has  been  obtained  by  an  exact 
solution  of  the  small  perturbation  transonic  equation,  using  guidelines  and  tools 
developed  at  the  Naval  Postgraduate  School.  This  unique  contour  allows  the  delay 
of  shock  formation  in  the  aft  portion,  hence  delaying  the  onset  of  wave  drag  which 
results  in  a  reduction  of  the  overall  transonic  pressure  drag  on  the  body.  The  drag 
coefficient  thus  computed  is  compared  with  another  axi-symmetric  body  with  the 
same  ellipsoid  forebody  but  a  simple  boat-tailed  conical  afterbody.  Computational 
Fluid  Dynamics  (CFD)  has  been  used  to  compute  the  viscous  flow  over  the  two 
bodies  using  a  Navier-Stokes  flow-solver.  Results  obtained  confirm  the  advantage 
of  the  special  shaped  afterbody  over  the  conical  afterbody  by  showing  the  delayed 
formation  of  shock  waves  at  the  aft  portion  in  transonic  flow,  consequently 
achieving  a  lower  maximum  drag  coefficient  of  approximately  5.5%.  These  results 
can  be  used  in  the  design  of  low  pressure-drag  surfaces  for  shapes  such  as 
missiles,  projectiles,  aircraft  external  ferry  tanks  and  aircraft  engine  nacelles  for 
improved  performance  within  the  transonic  flight  regime. 
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I .  INTRODUCTION 


The  study  of  the  aerodynamic  characteristics  of  a 
missile-like  object  from  launch  to  nominal  flight  speed  will 
typically  range  from  low  subsonic  to  supersonic  speeds. 
Entrenched  within  this  range  is  a  flight  regime  known  as  the 
transonic  flight  regime. 

The  study  of  transonic  aerodynamics  deserves  individual 
treatment  within  the  framework  of  gasdynamics.  A  typical 
transonic  flow  with  a  subsonic  free-stream  Mach  number 
contains  a  supersonic  zone,  bounded  by  the  sonic  line  and  a 
shock  wave,  through  which  the  flow  decelerates  back  to 
subsonic  flow.  Transonic  flows  occurs  typically  within  the 
free-stream  Mach  nximber  range  of  0.8  <  M«  <  1.2,  depending  on 
the  shape  of  the  object  [Ref.  1] . 

Transonic  aerodynamics  is  difficult  to  analyze  or  predict 
because  it  must  be  described  by  nonlinear  equations.  The 
transonic  flow  regime  is  also  complicated  by  the  formation  of 
shock  waves,  the  location  of  which  depend  on  the  comple:: 
interaction  of  free-stream  Mach  number  and  geometry  of  the 
body.  This  is  in  contrast  to  the  purely  subsonic  or  supersonic 
flow  regimes  where  an  adequate  prediction  of  aerodynamic 
characteristics  can  be  obtained  using  linear  theory. 
Consequently,  transonic  flows  are  very  sensitive  to  small 
perturbations  in  the  various  flow  param.eters  such  as  Mach 
number  and  this  makes  accurate  experimentation  difficult  [Ref. 
1]  . 

One  of  the  main  aerodynamic  forces  of  interest  is  drag. 
Objects  that  achieve  flight  are  designed  with  aerodynamic 
shapes.  The  prime  objective  of  such  designs  is  to  obtain  high 
lift  and  low  drag.  A  condition  associated  with  transonic  flow 
is  the  significant  increase  in  drag  [Ref.  2]  .  Efforts  to 
reduce  the  drag  coefficient  in  the  transonic  regime  must 
concentrate  on  reducing  the  pressure  drag  (inclusive  of  wave 


1 


drag)  contributions  to  total  drag.  During  wing  design,  one 
method  to  achieve  lower  drag  is  to  use  supercritical  airfoil 
sections  which  are  designed  to  delay  the  formation  of  a  shock 
wave  on  the  upper  surface  of  the  airfoil,  thereby  enabling  a 
higher  region  of  lift  generation  and  a  weaker  shock  to  occur 
further  downstream  of  the  airfoil  upper  surface  [Ref.  2] . 

A  traditional  approach  to  the  investigation  of 
aerodynamic  characteristics  is  based  on  wind  tunnel  test  data 
and  actual  flight  test  results.  Unfortunately,  both  wind 
tunnel  and  flight  testing  are  considerably  expensive  and  time 
consuming.  Furthermore,  wind  tunnel  testing  in  the  transonic 
regime  requires  special  wind  tunnels  with  adaptive  walls  in 
order  to  nullify  possible  reflected  shock-interaction  between 
the  tunnel  walls  and  the  test  object  within  the  test  section 
[Ref.  1].  As  mentioned  in  Reference  1,  it  is  also  difficult  to 
obtain  highly  reliable  data  from  wind  tunnel  experimentations 
in  the  transonic  flow  regime. 

In  contrast,  there  has  been  much  progress  recently  in  the 
field  of  Computational  Fluid  Dynamics  (CFD) .  The  judicious  use 
of  CFD,  in  combination  with  the  appropriate  modelling,  has 
enabled  aerodynamicists  to  predict  aerodynamic  characteristics 
within  a  faster  time  frame  and  at  relatively  lower  costs. 
Numerous  codes  are  now  available  that  run  as  both  the  Euler 
and  Navier-Stokes  flow-solvers  [Ref.  3]  .  Many  institutions  are 
now  using  CFD  as  a  "first-cut"  predictor  for  improvements  in 
aerodynamic  designs  or  as  a  tool  to  extrapolate  aerodynamic 
characteristics  from  known  design  performances.  With  the 
increased  experience  in  CFD,  such  as  better  turbulence 
modelling  and  more  efficient  and  robust  codes,  and  with  modern 
developments  in  high  speed  computing,  one  can  expect  to  see  an 
even  larger  proportion  of  aerodynamic  design  being 
accomplished  using  CFD. 
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II.  TRANSONIC  FLOW 


Transonic  flows  are  characterized  by  the  simultaneous 
presence  within  the  flow  field  of  both  subsonic  and  supersonic 
regions.  The  properties  of  transonic  flows  can  be  described 
using  the  equations  of  gas  dynamics,  namely,  the  equations  of 
state,  continuity,  momentum  and  energy.  If  the  object  of 
interest  is  a  slender  body,  we  can  use  small  perturbation 
theory  to  simplify  the  flow  over  the  body  [Ref.  2] .  In  fact, 
small  perturbation  theory  can  be  applied  throughout  the  entire 
Mach  number  range  from  low  subsonic  to  supersonic  speeds . 
Typically,  the  derivation  is  based  on  assumptions  that  the 
flow  is  steady,  irrotational  and  isentropic  with  no  energy 
transfer,  no  body  forces  and  no  shear  stresses  (i.e.,  inviscid 
flow) . 

A.  SMALL  PERTURBATION  THEORY 

For  steady,  isentropic  and  inviscid  flow,  the  governing 
equations  are 


Continuity  Equation 


d{pu)  ^d(pv)  ^  dipw) 
dx  dy  dz 


Expanding  and  re-arranging  equation  (1) 


p  3“  +  p  H- p  =  -  ( u #  +  v4£  +  ) 

^  dx  ^  dy  ^  dz  dx  dy  dz 


(1) 


(2) 
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ii)  The  Euler  momentum  equations  are  given  by 


du 

ox 

du 

dy 

du 

p 

dx 

dv 

dv 

.^dv_ 

dp 

dx 

dy 

dz 

p 

dy 

dw 

dw 

__i 

dp 

dx 

dy 

dz 

p 

dz 

Since  the  flow  is  isentropic,  we  can  express  the  speed  of 
sound  as 


Multiplying  the  x-momentum  equation  (Eqn  3)  by  u  yields 


7  du  du  du  u  dp  dp  in) 

dx  dy  dz  p  ox  op 

Substituting  Eqn  (6)  into  Eqn  (7)  yields 

u=4H.u„4y*u,.|M— (8) 

dx  dy  dz  p  ox 

Similarly,  by  multiplying  the  y-momentum  and  z-momentum 
equations  by  v  and  w  respectively  yields 


dv  7  dv  dv  V  dp 
uv-^+v2  41+w-=-=- - 


dx  dy  dz  p  dy 


dw  2  dw _  wdp  2 
dx  dy  dz  p  dz 
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Summing  up  Eqns  (8),  (9)  and  (10)  gives 

3v  .  __5  dv _ dv _ dw 


u 


dx  dy 


dw 


dw 

dz' 


-^(u4e*v4f»v|e) 

p  dx  dy  dz 

Substituting  Eqn  (2)  into  Eqn  (11)  yields 


(11) 


u 


dx  dy  dz  dx  dy  dz  dx  dy 


.dv 


.dw^,„,dw 


2  dw. 
dz 


'  dx  dy  dz 


By  rearranging  Eqn  (12) ,  we  obtain 

,  du  _  I  —  ) 

3x  dx'  dy  dy  dz  a^  dz 


,  uv  du  ^  uv  dv  V  _  /  uw  du  ^  uw  dw  ^ 
a^  dy  a^  dx  dz  a^  dx 


(12) 


-  ( J^.§X  +  —  —  )  =0 
a^  dz  a^  dy 

Since  the  flow  is  irrotational ,  we  have 

du  _  dv 
dy  dx 


(13) 


(14) 


du  _  dw 
dz  dx 


(15) 


dw_  dv 
dy  dz 


(16) 
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Substituting  Eqns  (14),  (15)  and  (16)  into  Eqn  (13 

rearranging  yields 

)  and 

dx  a  2  dy  a  ^  dz  a  ^ 

^  uv  du  ^  uw  dw  ^  vw  dv_f. 
a^dy  a2  dx  ^  a^  dz 

(17) 

For  a  slender  body  moving  in  the  x-direction, 
perturbation  theory  gives  the  velocity  field  as 

small 

u=U„*u'  }  v=v'  }  w=w' 

(18) 

Since 

u'>u'  ;  u>v'  ;  U'>w' 

(19) 

We  can  simplify  Eqn  (17)  to 

a  2  dx  dy  dz 

(20) 

iii)  The  flow  is  isentropic  and  the  energy  equation  can  be 

expressed  as 

=h + —  =  Cons  tan  t 

(21) 

We  have 

h  . 

“2  2 

(22) 

We  can  use  the  isentropic  relationships 

h=CpT  ;  C=-.  .  ;  a-=yRT  ;  M-  ^ 

^  -P  (y-l)  a 

(23) 
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Substituting  Eqns  (23)  into  Eqn  (22)  and  rearranging  yields 

a2=ai-I^  (V2-vf)  (24) 


v2-vf 


-)id 


Since  the  flow  is  3  dimensional,  we  have 

Substituting  in  Eqns  (18)  and  rearranging  yields 

V^-Ul^2Ujl'^u'^+v'^+w^ 

Substituting  Eqn  (27)  into  Eqn  (25)  yields 


^=1-1^  { 
^2  2 


2U^u'+u'^+v'^  +  w'^ 


Taking  the  reciprocal  of  Eqn  (2! 


2Ujj'+u'^+v'^+w'^ 


Using  the  Binomial  Expansion  on  Eqn  (29) 

■gj  O.  T. 


From  Eqn  (20),  the  first  term  can  be  expanded  as 


(l-^)-l-( 


Uf+2U„+U'^ 

=3  2  ^  ,,2  ^  2  ' 

a  u~a^ 


a^2U..u'^  ^ 

a  2 
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Substituting  Eqn  (30)  into  Eqn  (31)  yields 


U-  Uf 


u.  ul  ul 


Eqn  (32)  can  be  simplified  by  neglecting  the  higher  order 
terms  (H.O.T.)  to  give 

(l-ii!)  =1-Wf  [1  +  4^^1  +  -^A^)]  <33) 

a  2  z 

Substituting  Eqn  (33)  into  Eqn  (20)  yields 

=0  <34) 

Rearranging  Eqn  (34)  yields 

(1-Afi)  (35) 

dx  dy  dz  2  ax 

Rewriting  Eqn  (35)  in  terms  of  the  perturbation  velocities 


gives 


We  assume  that 


Expanding  Eqn  (36) ,  substituting  into  Eqn  (37)  and  neglecting 
2nd  order  perturbation  terms  yields 

(38) 

dx  ay  dz  2  ax 


Since  the  flow  is  irrotational,  we  have 


u'=- 


dy  dz 


(39) 


Substituting  Eqn  (39)  into  Eqn  (38)  gives  the  small 
perturbation  potential  equation  for  transonic  flows 


( 1  -id )  ^  ^  +  l^ld ) 


(40) 


B.  AN  EXACT  SOLUTION  FOR  AXI -SYMMETRIC  BODIES  IN  TRANSONIC 
FLOW 

Solutions  to  the  two-dimensional,  small  perturbation 
potential  equation  for  transonic  flow  have  been  accomplished 
by  Biblarz  [Ref.  4]  by  using  the  separation  of  variables 
approach.  The  solution  of  the  small  perturbation  transonic 
flow  equation  has  enabled  the  generation  of  shockless  boundary 
surfaces.  Al-Hashel  [Ref.  5]  reports  on  axi-syimnetric  boundary 
surfaces  for  M„=1.0,  1.10  and  1.20. 


9 


III.  DRAG 


Drag  can  be  defined  as  the  fluid  dynamic  force  acting 
against  an  object  placed  in  the  stream  of  airflow.  Basically, 
the  drag  force  can  be  divided  into  two  broad  categories.  The 
categories  are  viscous  or  skin  friction  drag  and  pressure 
drag . 

Viscosity  is  the  molecular  "resistance"  which  fluid 
particles  exhibit  against  displacement  in  relation  to  each 
other  and  with  respect  to  the  surface  of  solid  objects.  Such 
resistance  presents  itself  in  the  form  of  frictional  force, 
which  is  comparable  to  that  of  solid  surfaces  sliding  along 
each  other.  This  frictional  force  is  a  tangential  force  which 
occurs  when  a  fluid  (such  as  air)  flows  past  a  body.  This 
force  is  the  skin  friction  drag.  At  very  low  speeds  or  in 
flows  associated  with  fluids  of  high  viscosity  and  bodies  of 
very  small  dimensions  (i.e.,  at  very  low  Reynolds  number 
flow) ,  viscosity  is  the  predominant  property  determining  the 
drag  of  a  body  [Ref.  6].  At  higher  Reynolds  number  flows, 
other  forms  of  drag,  such  as  turbulent  skin  friction  drag, 
begins  to  have  significant  contribution. 

In  contrast  to  the  skin  friction  drag  (which  is  a 
tangential  force) ,  pressure  drag  results  from  the  distribution 
of  pressure  forces  normal  to  the  body  surface.  Within  the 
subsonic  range  of  speeds,  inviscid  fluid  dynamic  theory 
predicts  that  the  flow  will  close  in  behind  an  object,  without 
any  losses.  Positive  pressures  at  the  front  of  the  object  are 
counterbalanced  by  the  negative  pressure  of  equal  magnitude  on 
the  rear,  resulting  in  the  pressure  drag  being  zero.  This 
phenomenon  is  known  as  D'Alembert's  Paradox.  In  real  flows, 
however,  such  a  phenomenon  does  not  exist.  The  pressure  at  the 
rear  side  of  an  object  can  be  very  much  lower  than  predicted 
from  non-viscous  flow  conditions.  Viscosity  causes  the 
formation  of  a  boundary  layer  which  may  cause  flow  separation 


on  the  surface  of  an  object.  The  resultant  pressure 
differential  between  the  pressures  on  the  forward  side  and  the 
lower  pressures  on  the  rear  side  represents  the  pressure  drag. 
In  blunc  or  bluff  bodies,  this  type  of  drag  is  usually  much 
higher  than  the  skin  friction  drag  (Ref.  6].  For  streamlined 
or  aerodynamically  shaped  bodies,  the  reverse  is  true.  For 
flow  in  the  transonic  regime,  the  formation  of  shocks 
contributes  to  the  pressure  drag  due  to  the  pressure 
differential  between  the  forward  portion  of  an  object  and  the 
rear  portion  of  the  object.  Figure  (3.1)  [Ref.  7]  shows  the 
typical  drag  composition  of  an  axi-symmetric  projectile  shaped 
body.  Figure  (3.2)  [Ref.  8]  shows  the  variation  in  skin 
friction  drag  and  wave  drag  as  a  function  of  free-stream  Mach 
number  for  a  body  of  revolution  at  zero  incidence. 


Figure  3.1.  Components  of  Drag  [Ref.  7] 
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Figure  3.2.  Friction  and  Wave  Drag  Variation  [Ref.  8] 

One  well-known  method  of  reducing  the  drag  on  an  axi- 
symmetric  body  is  by  streamlining  the  afterbody,  i.e., 
reducing  its  base  diameter  gradually.  Figure  (3.3)  [Ref.  6] 
shows  the  reduction  in  zero  incidence  drag  coefficient  for 
boat-tailed  afterbodies.  Other  approaches  of  reducing  drag  on 
a  streamlined  body  includes  mechanical  methods  such  as  vortex 
generators  and  methods  for  controlling  the  boundary  layer  such 
as  the  boundary  layer  suction  method. 
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H  »  0.8  ;  Rtt  -  6  *  10’ 
Direction  of  Flow  - > 


Figure  3.3.  Drag  Coefficient  of  Boat-Tailed  Afterbodies  (Ref. 


Numerous  methods  are  available  to  calculate  the 
theoretical  values  of  drag  on  an  object.  However,  there  are  no 
totally  satisfactory  method  of  accurately  calculating  viscous 
drag,  particularly  transonic  drag.  This  is  due  in  large  to  the 
necessary  assumptions  and  simplifications  to  the  full  Navier- 
Stokes  equations  to  facilitate  mathematical  solutions  and  the 
need  to  use  turbulence  models  for  the  boundary  layer  effects 
in  high  Reynolds  Number  flows.  In  this  work,  both  the  skin 
friction  drag  and  the  pressure  drag  can  be  obtained  directly 
from  the  results  of  CFD  computations  [Ref.  3]. 
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IV.  COMPUTATIONAL  FLUID  DYNAMICS  (CFD) 


Rapid  advancement  in  the  speed  and  efficiency  of 
computing  and  better  modellings  has  led  to  the  emergence  of 
computational  fluid  dynamics  as  an  important  branch  of 
aerodynamics.  CFD  complements  the  experimental  and  theoretical 
branches  of  fluid  dynamics  by  providing  a  relatively  accurate 
and  cost-effective  means  of  simulating  real  flows.  Another 
significant  advantage  of  CFD  is  the  convenience  of  being  able 
to  vary  any  relevant  parameter  of  the  flow  and  also  the 
ability  to  "switch  off"  specific  terms  in  the  governing 
equations  so  as  to  assist  the  researcher  in  understanding 
their  contributions  to  the  resultant  flow. 

In  this  work,  CFD  was  used  to  compute  the  axi-symmetric 
flow  over  complete  body  models  which  are  composed  of  a 
forebody  (ellipsoid)  and  an  afterbody  (conical  or  special 
shape)  using  a  Navier-Stokes  flow-solver.  The  surface  profile 
of  the  special  shape  afterbody  was  reported  on  by  Al-Hashel 
[Ref.  5]  using  an  exact  solution  method  for  small 
perturbation.  The  conical  afterbody  has  a  base  diameter  ratio 
(d,./d,, )  of  0.75  and  a  conical  turning  angle  (P)  of  14  degrees. 
A  generic  forebody  in  the  shape  of  an  ellipsoid  was  patched 
onto  the  afterbodies  to  establish  the  CFD  working  models. 

A.  GRID  GENERATION 

A  computational  grid  system  is  a  necessary  part  of  any 
numerical  solution  based  on  a  finite-difference  or  finite- 
element  method.  The  selection  of  a  grid  system  is  based 
primarily  on  the  requirement  for  accuracy  in  the  final 
solution.  Secondary  considerations  are  the  effects  on 
computational  efficiency  of  the  solution  algorithm  using 
available  computer  architecture  and,  finally,  the  ease  of  grid 
generation. 
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The  grids  for  the  two  models  used  were  generated  by 
Priyono  [Ref.  9]  using  the  computer  programs  GRAPE  [Ref.  10] 
and  GRIDGEN2D  [Ref.  11] .  These  programs  were  used  to  generate 
two-dimensional  grids  which  were  subsequently  converted  into 
3-D  grids  utilizing  the  FORTRAN  codes  d2d3.f  and  rotategr.f 
[Appendix  B]  .  GRAPE  was  used  to  generate  the  grid  for  the 
model  with  the  special  shape  afterbody  (model  CNTRll) . 
GRIDGEN2D  was  used  to  generate  the  grid  for  the  model  with  the 
conical  afterbody  (model  CBXCN2),  since  GRAPE  cannot 
accomodate  the  non-smooth  points  associated  with  the  conical 
afterbody.  The  grids  for  the  two  models  are  shown  at  Figures 
(1)  and  (2)  in  Appendix  A.  In  this  work,  the  correspondence 
between  the  two  body  profiles  was  verified  prior  to  further 
analysis.  This  correspondence  is  shown  in  Figure  (3)  in 
Appendix  A. 

Since  the  analysis  was  for  the  complete  body,  an  0-type 
grid  was  used.  Previous  work  by  Priyono  demonstrated  the  need 
to  use  a  fine  grid  due  to  the  significant  propagation  of  the 
bow  shock  above  and  below  the  body  in  transonic  flow, 
especially  in  the  vicinity  of  freestream  Mach  number  of  unity. 
It  was  necessary  for  the  bow  shock  to  diminish  before  it 
reached  the  grid  boundary,  thereby  ensuring  that  the  shock  is 
not  reflected  by  the  grid  boundary  and  contaminated  the 
solution.  The  outer  boundary  of  the  fine  grid  was  5  times  the 
body's  length.  The  grids  were  generated  only  for  half  of  the 
complete  body  taking  advantage  of  symmetry.  The  grid  size  used 
IS  77x15x120. 

B.  FLOW- SOLVER 

The  OVERFLOW  program  was  developed  by  NASA  Ames  Research 
Center.  The  version  used  is  Version  l.Sag  released  on  30  April 
1993  [Ref.  3].  The  program  uses  either  the  3-D  Euler  or 
Navier-Stokes  flow-solvers  for  inviscid  or  viscous  flow 
computations  respectively.  The  selection  of  Euler  or  Navier- 
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stokes  flow-solver  is  accomplished  by  setting  the  parameter 
VISINP  (viscosity  input)  to  True  or  False  in  the  input  file. 

Before  applying  the  flow-solver  code  (OVERFLOW),  the 
formatted  3-D  grid  file  (grid. for)  must  be  converted  into  an 
unformatted  input  grid  file  (grid. in) .  This  is  accomplished 
using  the  FORTRAN  code  readx.f  [Appendix  B]  .  An  input  file 
(overflow. in)  containing  all  the  input  parameters  must  then  be 
written  for  running  OVERFLOW.  The  input  parameters  consist  of 
the  number  of  iterations,  timesteps,  calculation  methods, 
smoothing  method,  types  of  in-flow  and  out-flow,  boundary 
conditions  and  types  of  turbulence  models.  The  value  of  angle 
of  attack  (ALPHA)  depends  on  the  orientation  of  the  grid  in 
the  coordinate  system.  In  our  case,  ALPHA  was  set  at  180“ 
(i.e.,  flow  comes  from  the  x-positive  to  the  x-negative 
direction).  A  typical  input  file  (overflow . in)  is  provided  in 

Appendix  B. 

In  his  work,  Priyono  investigated  the  flow  past 
afterbodies  only  using  both  the  Euler  and  Navier-Stokes  flow- 
solvers.  Priyono  also  investigated  the  flow  past  complete 
bodies  but  only  with  the  Euler  flow-solver.  For  the  present 
work  which  investigates  the  flow  past  complete  bodies  using 
the  Navier-Stokes  flow-solver,  the  effects  of  viscosity  have 
to  be  “switched  on”.  This  is  accomplished  by  setting  VISCJ, 
VISCK  and  VISCL  parameters  to  TRUE.  This  in  effect  includes 
viscous  terms  into  the  J  (along  the  body  in  the  x-direction) , 
K  (around  the  body  or  x-axis)  and  L  (outwards  from  the  body 
towards  the  outer  grid  boundary)  computational  directions.  The 
number  of  turbulent  wall  regions  NTURB  is  set  to  2.  This 
corresponds  to  the  turbulence  models  ITTYP  of  1  and  11,  i.e., 
the  Baldwin-Lomax  boundary  layer  model  (with  variable  Degani- 
Schiff  cutoff)  for  the  body  in  general  and  the  Baldwin-Lomax 
shear  layer  model  (with  variable  C,,,)  to  model  the  wake  behind 
the  body.  The  basis  of  selecting  the  turbulence  models  comes 
from  the  OVERFLOW  manual  which  recommends  the  use  of  the 
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Baldwin -Lomax  model  for  flows  with  mild  separation. 

The  OVERFLOW  program  generates  output  files  such  as 
overflow. out,  fomo.out,  rpmin.out,  resid.out  and  q.save.  The 
file  overflow. out  provides  information  on  the  input  parameters 
which  are  used  in  the  computational  run.  The  file  fomo.out 
provides  the  force  and  moment  history,  which  in  this  case, 
includes  coefficients  of  drag  due  to  pressure  (which  was  also 
computed  when  using  the  Euler  flow  solver)  and  viscosity 
(i.e.,  skin  friction  drag).  The  file  rpmin.out  provides  the 
minimum  density/pressure/ratio  of  specific  heats  history.  The 
file  resid.out  provides  flow  solver  residual  history  for  the 
five  primary  variables  (i.e.,  p,  pu,  pv,  pw  and  e,  where  e  is 
the  total  energy) .  The  variable  that  is  used  to  show 
convergence  of  the  solution  is  the  density  residual.  In 
general,  the  convergence  criterion  can  be  defined  as  the 
reduction  in  the  density  residuals  by  two  orders  of  magnitude 
and  the  trend  towards  further  reduction.  The  file  q.save  is 
the  solution  file  in  PLOT3D  format  [Ref.  12]  when  OVERFLOW  has 
run  satisfactorily.  If  the  flow-solver  detects  a  case  of 
negative  density  or  pressure  during  the  computation  process, 
an  output  file  called  q.bomb  will  be  generated  instead  of  the 
solution  output  file  q.save.  The  file  q.bomb  thus  represents 
the  solution  prior  to  the  error  being  detected. 

The  output  from  resid.out  is  a  text  file  and  the  result 
can  be  plotted  using  a  plotting  software  such  as  GNUPLOT 
[Ref.  13].  The  output  from  fomo.out  is  also  a  text  file  and 
the  values  of  coefficient  of  drag  due  to  pressure  and 
viscosity  can  be  read  off  directly  at  the  end  of  the  file.  The 
coefficient  of  drag  can  also  be  plotted  using  GNUPLOT  to  shov; 
convergence  to  a  particular  value.  The  output  from  q.save  is 
an  unformated  file  which  needs  to  be  converted  into  a 
formatted  file  (q.form)  using  the  FORTRAN  code  readq.f 
[Appendix  B]  .  The  formatted  solution  file  q.form  can  be 
plotted  using  either  PL0T3D  or  FAST  [Ref.  14].  Both  software 
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packages  are  capable  of  generating  plots  showing  the  Mach 
number  distribution,  shock  locations  due  to  pressure 
difference,  velocity  vector  distribution  and  many  other 
options  of  flow  visualization  over  the  bodies.  In  this  work, 
PL0T3D  was  used  for  all  the  flow  visualization  plots. 

C .  CPD  SOLUTIONS 

Initial  turbulent  flow  computations  were  accomplished 
with  the  setting  of  ITDIR  of  3  and  1,  corresponding  to  the  L 
and  J  coordinate  directions  away  from  the  wall  or  shear  layer 
for  the  Baldwin-Lomax  boundary  layer  model  and  the  Baldwin- 
Lomax  shear  layer  model  respectively.  Solutions  obtained  from 
the  computations  displayed  an  asymmetry  at  the  wake  of  the 
bodies.  Logically,  there  should  be  no  asymmetry  since  the 
computations  were  for  zero  incidence  flow. 

Re-examination  of  the  input  for  ITDIR  strongly  indicated 
that  the  parameter  should  be  set  to  3  and  3,  corresponding  to 
the  coordinate  direction  of  L  for  both  the  turbulence  models. 
Subsequent  computations  with  this  setting  produced  symmetry  at 
the  wake.  The  output  also  displayed  a  clear  momentum  deficit 
within  the  wake  region  aft  of  the  bodies. 

Having  resolved  the  parameters  for  the  turbulence 
modelling,  the  next  step  was  to  ensure  effective  convergence 
of  the  solutions.  A  slow  start  was  used  with  DT  =  0.01  and 
CFLMIN  =  3.0  for  all  the  runs.  After  2000  iterations,  it  was 
observed  that  the  density  residuals  were  only  converging  to 
appi'oximat  ely  1.5  orders  of  magnitude.  Since  the  viscous 
solution  includes  skin  friction  drag,  the  residual  would  have 
to  show  at  least  2  orders  of  magnitude  of  convergence.  In 
order  to  increase  the  convergence  rate,  the  runs  were 
initiated  for  either  500  or  1000  iterations  with  a  CFLMIN  =  3, 
depending  on  the  free-stream  Mach  number.  Subsequently, 
restarts  were  used  and  the  values  of  CFLMIN  or  DT  were  varied. 
Typical  variations  of  CFLMIN  and  DT  are  tabulated  at  Table  (1) 
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in  Appendix  C.  Figures  (4)  and  (5)  in  Appendix  A  show  typical 
plots  for  density  residuals.  The  fluctuations  in  the  residuals 
are  caused  by  the  effects  of  turbulence  modelling. 

The  critical  free-stream  Mach  number  for  both  bodies  was 
computed  to  be  approximately  past  0.8.  This  is  in  contrast  to 
the  Euler  results  reported  by  Priyono  which  show  the  inviscid 
critical  Mach  number  as  approximately  0.7  and  that  sonic 
conditions  were  reached  in  the  regions  where  the  forebody  and 
afterbodies  mate  with  the  mid-sections .  Figures  (6)  and  (7)  in 
Appendix  A  shows  the  Mach  contours  for  CNTRll  and  CBXCN2 
respectively  at  M„=0.8  and  the  highest  local  Mach  number 
displayed  is  0.9,  indicating  that  the  critical  Mach  number  has 
not  been  reached.  Figures  (8)  and  (9)  in  Appendix  A  shows  the 
Mach  contours  and  shock  locations  for  CNTRll  at  M»=0.95 
respectively.  Comparing  Figure  (9)  with  Figure  (10)  in 
Appendix  A,  which  shows  the  shock  locations  for  CBXCN2  at 
M„=0.95,  note  that  the  rear  shock  is  further  back  at  the  aft 
portion  of  CNTRll.  The  viscous  critical  Mach  number  was  higher 
due  to  the  presence  of  the  boundary  layer  which  lessened  the 
acceleration  of  the  airflow  over  the  forebody  and  also  had  a 
"rounding"  effect  on  the  afterbodies  which  translates  into  a 
slower  acceleration  around  the  aft  regions. 

Since  this  work  uses  the  Navier-Stokes  flow— solver,  it  is 
important  to  verify  that  the  resultant  flow  is  indeed  a 
viscous  solution.  The  plot  at  Figure  (11)  in  Appendix  A  shows 
the  velocity  vector  distribution  over  CNTRll  at  a  free-stream 
Mach  number  of  1.2.  It  is  important  to  note  the  distinct 
boundary  layer  in  the  plot.  It  was  also  necessary  to  validate 
the  no-slip  condition  on  the  surface  of  the  bodies.  The 
appropriate  input  parameters  were  used  in  PLOT3D  to  show  only 
the  Mach  contours  on  the  surface  of  the  bodies.  The  results 
show  a  Mach  contour  of  zero  on  the  surface  of  the  bodies, 
thereby  confirming  the  no-slip  condition. 
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The  plots  showing  the  shock  locations  based  on  pressure 
gradient  complement  the  plots  of  Mach  number  distribution.  The 
plots  at  Figures  (12)  and  (13)  in  Appendix  A  for  CNTRll  and 
CBXCN2  respectively  at  M,=1.05  show  that  a  bow  shock  has 
formed  once  the  free-stream  Mach  number  goes  supersonic.  It  is 
interesting  to  note  that  the  plot  at  Figure  (10)  in  Appendix 
A  for  free-stream  Mach  number  of  0.95  corresponds  to  a 
schlieren  picture  at  Figure  (14)  in  Appendix  A  of  a  similar 
flow  over  an  axi-symmetric  body  at  free-stream  Mach  number  of 
0.95  at  zero  incidence  [Ref.  15]. 
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V.  RESULTS,  CONCLUSIONS  AND  RBCOMMENDATIONS 


A.  RESULTS 

The  values  of  drag  coefficient  are  obtained  directly  from 
Che  OVERFLOW  output  files  fomo.out.  The  output  from  these 
files  iS  plotted  to  show  that  the  values  of  pressure  drag 
coefficient  and  skin-friction  drag  coefficient  are  converged. 
A  typical  plot  is  shown  at  Figure  (5.1)  .  The  viscous  drag 
coefficients  of  the  two  bodies  at  the  various  Mach  numbers  are 
compared  in  Figure  (5.2) .  The  corresponding  inviscid  drag 
coefficients  as  computed  by  Priyono  are  indicated  in  Figure 
(5.3). 


CBXCN2  (NAVIER  STOKES  SOLUTION) 


Figure  5.1.  Convergence  of  Drag  Coefficient 
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Fr«»  Stream  Mach  Number 

Figure  5.2.  Drag  Coefficient  For  Navier- Stokes  Solutions 


Cd  v8  Mach  Number  (Euler  Solution) 


As  expected,  the  viscous  drag  coefficients  are  higher 
than  the  inviscid  drag  coefficients.  These  results  show  that 
skin  friction  drag  coefficients  are  significant  for  the 
viscous  solutions.  The  viscous  solutions  show  a  maximum  drag 
coefficient  difference  of  approximately  5.5%  between  the  two 
bodies.  In  comparison,  the  Priyono  [Ref.  9]  reported  a  maximum 
drag  coefficient  difference  of  approximately  15%  between  the 
two  bodies  for  the  Euler  solutions.  The  likely  difference 
between  these  is  the  smoothing  of  the  afterbody  by  the 
boundary  layers . 

The  pressure  drag  coefficient  remains  relatively  constant 
from  the  sub-sonic  free-stream  Mach  number  regime  up  until  the 
critical  Mach  number.  Subsequently,  the  pressure  drag 
increases  due  to  the  addition  of  wave  drag  resulting  from  the 
formation  of  shock  waves.  The  skin  friction  drag  coefficient, 
on  the  other  hand,  decreases  with  increasing  free-stream  Mach 
number.  The  drag  coefficients  as  shown  in  Figure  (5.2)  is  in 
good  agreement  with  typical  values  as  shown  in  Figures  (3.1) 
and  (3.2)  .  Preliminary  results  for  free-stream  Mach  numbers  of 
1.3  and  1.5  show  that  the  drag  coefficient  decreases  from  the 
peak  value  shown  at  Figure  (5.2). 

The  results  show  that  the  special  shape  afterbody  did 
indeed  provide  a  lower  drag  coefficient  as  compared  to  the 
conical  afterbody.  The  plots  from  shock  location  for  free- 
stream  Mach  number  of  0.95  clearly  shows  the  position  of  the 
rear  shock  wave  as  being  further  aft  of  the  special  shape 
afterbody  when  compared  with  the  rear  shock  wave  location  for 
CBXCN2 . 

B.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  OVERFLOW  code  has  proven  to  be  a  useful  CFD  tool  for 
this  work.  In  hindsight,  it  would  have  been  preferable  to 
generate  the  grids  for  both  the  bodies  using  one  grid 
generation  method  instead  of  using  GR.A.PE  for  CNTRll  and 
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GRIDGEN2D  for  CBXCN2 .  The  advantage  would  have  been  complete 
commonality  between  the  grids  for  the  two  bodies  with  the  only 
variations  being  at  the  aft  portions.  As  mentioned  earlier, 
although  the  profile  of  the  two  bodies  were  verified,  the 
positioning  of  the  grid  points  on  the  two  bodies  (especially 
on  the  forebody)  were  different.  This  difference  may  affect 
the  accuracy  of  the  solutions. 

The  investigation  of  the  two  axi-symmetric  bodies  using 
CFD  has  shown  an  advantage  of  the  special  shaped  afterbody  in 
the  form  of  lower  transonic  viscous  drag  coefficient  over  the 
conical  shaped  afterbody  by  a  maximum  difference  of 
approximately  5.5%.  This  advantage  is  over  and  above  the 
benefits  that  one  can  achieve  by  boat-tailing  the  afterbody  as 
shown  in  Figure  (3.3)  .  These  results  can  be  used  to  design  low 
pressure  drag  surfaces  for  shapes  such  as  missiles, 
projectiles,  aircraft  external  ferry  tanlcs  and  aircraft  engine 
nacelles  for  improved  performance  within  the  transonic  flight 
regime . 

Finally,  this  worJc  may  be  continued  with  experimental 
validation  of  the  advantage  of  the  special  shaped  afterbody 
over  the  simple  conical  afterbody  using  an  appropriate 
transonic  wind-tunnel.  Further  CFD  worlc  could  also  be  carried 
out  to  investigate  the  effects  of  angle-of -attack  on  the  two 
bodies.  Such  efforts  would  contribute  useful  experimental 
data  on  the  transonic  flow  regime,  the  likes  of  which  is 
certainly  lacking  currently. 
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APPENDIX  B.  FORTRAN  PROGRAMS  AND  OVERFLOW  INPUT  FILE 
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*  D2D3.F 

*  This  FORTRAN  program  will  pretJate  a  formatted  2-D  grid 

*  data  file  for  processing  by  ROTATEGR.F 


real  X(77, 120) , Y (77 , 120) , Z (77 , 120) 
open(unit=12,  f lle= 'grid2d.  fot ' ,  status= 'unlcnown' ) 
open (unit=14, f ile= 'cone3d. for' , status= 'new' ) 

READ (12,*)  IDIM,JDIM 

READ (12, *)  ( (X(I, J) , 1=1, IDIM) , J=1, JDIM) , 

+  ( (Y(I, J) ,1=1, IDIM) , J=l, JDIM) 

Z (I, J)=0 .0 

WRITE (14,*)  IDIM,JDIM,1 

WRITE ( 14, *)  ( (X(I, J) , 1=1, IDIM) , J=l, JDIM) , 

+  ( (0.0, 1  =  1, IDIM) ,J=1, JDIM)  , 

+  ( (Y(I,J) , 1=1, IDIM) ,J=1, JDIM) 

.STOP 

END 


*  ROTATEGR . F  * 

*  This  FORTRAN  program  talces  a  formatted  2-D  grid  and  * 

*  generates  an  Axi-Symmetric  an  unformatted  3-D  grid  * 

*  Converts  grid2d.for  to  cone. for  * 


dimension  x(77,120),  y(77,120),  z(77,120) 

dimension  xx { 77 , 15 , 120 ) ,  yy ( 77 , 15 , 120 ) ,  zz ( 77 , 15 , 120 ) 
character*30  fni 
character*30  fno 

print  *,  'Input  2-D  grid  filename' 
read  (*,21)  fni 
rewind  2 

open  (2, file=fni, f orm= ' formatted ' ) 
read  (2,*)  il,jl,kl 

read  (2,*)  ((  x ( i , j ) , i  =  l , il )  ,  j  =  l ,  j 1 )  , 

+  ((  y(i,j),i=l,il),j=l,jl), 

+  {(  z  (  i , j ) , i  =  l , il )  ,  3  =  1 ,  jl ) 

close (2 ) 

pi  =  4 . *atan ( 1 . ) 

print  *,  'No  of  Planes  in  j  dir  ?  ' 
read  ( 5 , *  )  jm 

dth  =  (180.  /  (jm-D)  *  (pi/180.) 

do  11  i  =  1, il 
do  11  j  =  l,jl 
k-j 

xx(i,2,)c)  =  x(i,k) 
y:>-  (  i  ,  2  ,  k )  =  0  . 
zz(i,2,k)  =  z(i,k) 

.11  continue  4  6 


im=il 

kin=jl 

do  20  j=3,jin+l 

do  20  i=l,iin 

do  20  k=l,km 

XX (i, j , k)  =  x(i, k) 

th  =  (j-2)*dth 

yy(i,j,k)  =  sin ( th) *z (i , k) 

2z{i,j,k)  =  cos ( th) *z { i , k) 

20  continue 

do  30  i=l,im 
do  30  k=l,km 
3  =  1 

xx{i,j,k)  =  xx(i,j+2,k) 
yy(i,j,k)  =-yy ( i , j +2 , k) 

2z(i,j,k)  =  2z(i,j+2,k) 
j  =  jin+2 

xx(i,j,k)  =  xx(i,j-2,k) 
yy ( i , j , k )  =  -yy ( i , j -2 , k ) 
zz(i,j,k)  =  zz(i,i-2,k) 

30  continue 

jm  =  jm+2 

print  *,  'Output  3-D  grid  filename  =' 
read  (5,21)  fno 
rewind  3 

open  (  3,  f ile-f no , form= ' unformat  ted ' ) 
write  (3)  im,im,km 

write  (3)  (((  xx(i,j,k),  i=l,im),  j=l , jm) , k=l , km) 

+  (((  yy(i,j,k),  i  =  l,im),  j  =  l ,  jm)  ,  k=l ,  )<m) 

+  (((  zz(i,j,k),  i  =  l,im),  i  =  l ,  jm)  ,  k=l ,  )an) 

close  (3) 
stop 
end 
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**************************************************************** 

*  READX.F  ,  * 

*  This  FORTRAN  program  will  read  formatted  data  files  and  * 

*  convert  them  to  unformatted  binary  data  filfs  for  input  * 

*  into  OVERFLOW  -  Converts  grid. for  to  grid. in  * 

**************************************************************** 


real  x(77,15,120) , y (77, 15, 120) , z (77, 15, 120) 
open (unit=12, file=' grid. for' , status= 'unknown' ) 

open (unit=14, file= 'grid. in' ,status= 'new' , form= 'unformatted' ) 
read (12,*)  IDIM, JDIM, KDIM 

read(12,*)  ( ( (X(I,J,K) ,1=1, IDIM) ,J=1,JDIM) ,K=1,KDIM) , 

(  (  (Yd,  J,K)  ,1=1, IDIM)  ,  J=l,  JDIM)  ,K=1,KDIM)  , 

( ( (Z(I, J,K) ,1=1, IDIM) , J=l, JDIM) ,K=1,KDIM) 
write (14)  IDIM, JDIM, KDIM 

write(14)  ( ( (X(I,J,K) ,1=1, IDIM) ,J=1, JDIM) ,K=1,KDIM) , 

( ( (Y(I, J,K) , 1=1, IDIM) , J=l, JDIM) ,K=1,KDIM) , 
(((Z{I,J,K), 1=1, IDIM) , J=l, JDIM) ,K=1,KDIM) 

stop 

end 


**************************************************************** 


*  READQ . F 

*  This  FORTRAN  program  will  read  unformatted  binary  data  * 

*  files  and  convert  them  to  formatted  data  files  for  input 

*  into  PLOT3D  -  Converts  q.save  to  q.form  * 

*-,r************************************************************** 


dimension  q(77,15,120,5) 

open (uni t=l, f ile= ' q . save ' , status= 'unknown' , form= 'unformatted' ) 
open(un  Lt=20, f ile= ' q. form' ,  status='new' , f orm= ' formatted' ) 
read ( 1 '  ni ,  nj ,  nk 
read ( i ) f smach, alpha , re , time 

read(l) ( ( ( (q(i, j , k,nx) , i=l,ni) , j=l,nj ) , k=l,nk) ,nx=l, 5) 

write (20,*)  ni,  nj,  nk 

write (20, * ) f smach, alpha, re, time 

write (20, * ) ( ( ( (q ( i , j , k, nx) , i=l , ni ) , j =1 , nj ) , k=l , nk) ,nx=l, 5) 

print*,  ni,  n j ,  nk 

stop 

end 
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INPUT  FILE  FOR  OVERFLOW  ( OVERFLOW . IN ) 
★************************************ 

$GLOBAL 

CHIMRA  =  .F.,  NSTEPS  =  1000,  RESTRT  =  .F.,  NSAVE  =  100, 
NQT  =  0 , 

$END 

$FLOINP 

ALPHA  =  180.0,  FSMACH  =  0.95,  REY  =  6.00E6,  TINF  =  520.0, 
$END 
$VARGAM 
$END 
$GRDNAM 

NAME  =  'Half  conplete  body,  CNTRll  ,  77x15x120  Grid', 

$END 

$NITERS 

$END 

$METPRM 

IRHS  =  0,  ILHS  =  2,  IDISS  =  2, 

$END 

$TIMACU 

DT  =0.01,  ITIME  =  1  TFOSO  =  1.00,  CFLMIN  =  3.0, 

$END 

$SMOACU 

ISPECJ  =  2,  DIS2J  =  2.00,  DIS4J  =  0.02, 

ISPECK  =  2,  DIS2K  =  2.00,  DIS4K  =  0.02, 

ISPECL  =  2,  DIS2L  =  2.00,  DIS4L  =  0.02, 

SMOO  =  1.00, 

EPSE  =  0.35, 

$END 

$VISINP 

VISCJ  =  .T.,  VISCK  =  .T.,  VISCL  =  .T., 

NTURB  =  2 , 

ITTYP  =  1,  11, 

ITDIR  =  3,  3, 

JTLS  =  -1,  1, 

JTLE  =  1,  -1, 

KTI.S  1,  1, 

KTLE  =  15,  15, 

LTLE  =  1 ,  1 , 

LTLE  =  -1,  -1, 

TLPAR1=  0.9,  1, 

$END 

$BCINP 


NBC 

= 

6, 

IBTYP 

= 

15, 

15, 

12, 

12, 

5, 

32, 

IBDIR 

= 

1, 

-1, 

2, 

-2, 

3, 

-3, 

JBCS 

= 

1, 

77, 

1, 

1, 

1, 

1, 

JBCE 

= 

1, 

77, 

77, 

77, 

77, 

77, 

KBCS 

= 

1, 

1, 

1, 

15, 

1, 

1, 

KBCE 

= 

15, 

15, 

1, 

15, 

15, 

15, 

LBCS 

=: 

1, 

1, 

1, 

1, 

1, 

120, 

LBCE 

$END 

;SCEINP 

120, 

120, 

120, 

120, 

1, 

120, 

$END 
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APPENDIX  C.  TABLE  OF  OVERFLOW  INPUTS 
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INITIAL  DISTRIBUTION  LIST 


No.  Copies 

1.  Defense  Technical  Information  2 

Cameron  Station 

Alexandria,  Virginia  22304-6145 

2.  Library,  Code  52  ^ 

Naval  Postgraduate  School 

Monterey,  California  93943-5101 

3.  Chairman  ^ 

Department  of  Aeronautics,  Code  AA/Co 

Naval  Postgraduate  School 
Monterey,  California  93943-5101 

4.  Prof.  0.  Biblarz,  Code  AA/Bi  3 

Naval  Postgraduate  School 

Monterey,  California  93943-5101 

5.  Prof.  G.  Hobson,  Code  AA/Hg  1 

Naval  Postgraduate  School 

Monterey,  California  93943-5101 

6.  Major  Fan,  Yue  Sang  2 

Air  Logistics  Department 

HQ  Republic  of  Singapore  Airforce 
MINDEF  Building,  Gombak  Drive 
SINGAPORE  2366 
Republic  of  Singapore 


53 


